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Long-Term Bioclimatic Modelling the Distribution of the Fire-Bellied Toad, Bombina bombina 
(Anura, Bombinatoridae), under the Influence of Global Climate Change. Tytar, V., Nekrasova, O., 
Pupina, A., Pupins, M., Oskyrko, O. — The article describes the potential distribution area of B. bombina 
and figure out the significant climatic factors of the species at a home range scale. This species is listed 
on Appendix II of the Bern Convention and on Annexes II and IV of the EU Natural Habitats Directive. 
It is protected by national legislation in many countries, occurs in many protected areas, and is listed in 
many national and sub-national Red Data books and lists. We collected the occurrence records and a 
set of climatic variables including 19 factors from 10’ resolution historical (summarizing annual trends, 
seasonality and extreme conditions during 1961-1990) and projected data (2050) available at the CliMond 
database. As a result, under climate predictions for 2050 there may be a substantial north and north-west 
shift of optimal habitat. Under such a scenario B. bombina populations may suffer mostly in the east and 
south of Ukraine. Under the modelled scenario the species representation in protected areas throughout 
the home range should be considered, but especially in Ukraine. 
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Introduction 
Recent reports suggest that as much as a third of all known amphibians are in decline, many inhabiting 


areas far from obvious human disturbances (Stuart, 2008). Global warming and climate change have been 
implicated as forces likely to drive amphibian declines by significantly changing a habitat through time (e. g., 
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Pounds, 2001; Carey, Alexander, 2003; Raffel et al., 2013). Meta-analysis comparing responses of a wide range 
of different taxonomic groups to climate change across several biogeographical regions already indicate shifts 
in the distribution patterns of many plants and animals (e. g., Parmesan, Yohe, 2003). The range shift parallels 
a 0.8 °C warming over Europe during the last century, which has shifted the climatic isotherms northwards by 
an average of 120 km (Beniston, 1998). The Intergovernmental Panel on Climate Change (Solomon, 2007 ;Cli- 
mate Change, 2014) projections forecast changes in the global climate during the 21st century even larger than 
those observed during the 20th century. Further temperature warming and decreasing of water availability will 
produce dramatic consequences for amphibians in terms of large extinctions and/or range shifts (Araújo et 
al., 2006). In the context of future climate change, range shifts are a key response, and can affect species repre- 
sentation in protected areas. Ectotherms are more likely to track their climate space compared to endotherms 
(Aragon et al., 2010) and major shifts in herpetofaunal assemblages caused by climate change are predicted 
worldwide (Lawler et al., 2010). Northward range shift would thus be reflected in either a net extinction at the 
southern boundary or a net colonization at the northern boundary. Because habitat suitability is shifting with 
climate change, such shifts could lead to population declines and high extinction risk for species that are unable 
to move to new appropriate conditions because of the patchy landscape and insufficient dispersal corridors 
(Araújo et al., 2006). 

The complex interaction between species distribution and climate change is yet insufficiently understood, 
but the increasing availability of information on the variation of environmental parameters in geographic 
space, species distribution data, and computation capacities during the last decade now allow large scale assess- 
ments of such relationships (Kozak et al., 2008). Relationships can be assessed by calculating ‘environmental’ 
or ‘ecological’ niches and their subsequent projection into geographic space (Guisan, Zimmerman, 2000). But 
testing the hypothesis that global warming has already caused amphibian declines is challenging because many 
factors could act synergistically with climate in complex ways (Pounds, 2001). 

In this study, we aimed at projecting the availability of suitable habitat for an endangered amphibian 
species, the Fire-bellied toad Bombina bombina (Linnaeus, 1761). The species is distributed in lowlands over 
a wide continental area: it is widespread in Poland, Lithuania, Belarus, Western Russian Federation, Ukraine, 
Bulgaria, Rumania and Hungary. The western border reaches the lowlands of Austria, Czech Republic and 
East-Germany. The northern border of distribution follows approximately the 56° degree of latitude. The most 
Northern European populations are in Southern Sweden, East-Denmark, Northern Germany (Schleswig-Hol- 
stein), Belarus and Southern Latvia (Drobenkov et al., 2005; Kuzmin et al., 2008). In Sweden during the 1960s, 
it became extinct. 

In western and northern Europe the species is threatened by the loss of habitat. Recent declines in north- 
western Europe might also be related to climate change (Agasyan et al., 2009). Even small climatic changes 
could have profound effects on the abundance of the species in this area. In Latvia, for instance, the cold climate, 
snowless winters, cold short springs and dry hot summers are specified as negative factors based on the risks for 
B. bombina populations (Pupina, Pupins, 2016). 


Material and methods 


A substantial amount of species records of B. bombina are available through the Global Biodiversity In- 
formation Facility (GBIF, www.gbif.org) and HerpNet databases (www.herpnet.org). In addition, species re- 
cords can be obtained from own field trips (Pupina, Pupins, 2015), museum collections or literature (Sillero 
et al., 2014). In total, 2303 georeferenced occurrence records of B. bombina finds were analyzed. Of these, 163 
are from Latvia (Pupina, Pupins, 2015) and 866 from Ukraine (80 % personal data of O.Nekrasova: Tytar, 
Nekrasova, 2016; Tytar et al., 2018 a, b and from the literature: Shaitan, 1999, etc.) that cover the southern and 
northern boundaries of the home range of B. bombina. 

In recent years, species distribution modeling (SDM) has been widely used to estimate ecological require- 
ment of particular species and to characterize and map the spatial distribution of habitat occupied by species 
at a landscape scale (Elith, Leathwick, 2009; Franklin, 2009; Peterson et al., 2011; Phillips, Elith, 2013). The 
principle of SDM is related to Hutchinson super-volume theory, and emphasizes the ecological requirement 
of species, especially the abiotic factors controlling species distribution (Guisan, Thuiller, 2005; Li et al., 2013). 
In general, many climatic factors were used as predicted variables when simulating species potential distribu- 
tion at large scale with coarse resolution. According to the predicted map, we can depict the climatic niche and 
response curves of species. Currently, there are many algorithms in SDM technique, but Elith et al. (2006) have 
shown that maximum entropy model (MaxEnt) is one of the best method among an array of algorithms. 

The aim of this work was to simulate the potential distribution area of B. bombina and figure out the 
significant climatic factors of the species at a home range scale. Firstly, we collected the occurrence records. 
Secondly, we collected a set of climatic variables including 19 factors from 10° resolution historical (sum- 
marizing annual trends, seasonality and extreme conditions during 1961-1990) and projected data (2050) 
available at the CliMond database (Kriticos et al., 2012) (table 1). These variables represent annual trends 
(e. g., mean annual temperature, annual precipitation) seasonality (e. g., annual range in temperature and 
precipitation) and extreme or limiting environmental factors (e.g., temperature of the coldest and warmest 
month, and precipitation of the wet and dry quarters). Finally, we used the occurrence records and cli- 
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Table 1. Description of bioclimatic variables 


on Variable 
Bio01 Annual mean temperature ,°C 
Bio02 Mean diurnal temperature range (mean; period max-min), °C 
Bio03 Isothermality, Bio02+Bio07 
Bio04 Temperature seasonality, C of V 
Bio05 Max temperature of warmest week, °C 
Bio06 Min temperature of coldest week, °C 
Bio07 Temperature annual range (Bio0-Bio06), °C 
Bio08 Mean temperature of wettest quarter, °C 
Bio09 Mean temperature of driest quarter, °C 
Biol0 Mean temperature of warmest quarter, °C 
Bioll Mean temperature of coldest quarter, °C 
Bio12 Annual precipitation, mm 
Bio13 Precipitation of wettest week, mm 
Biol4 Precipitation of driest week, mm 
Bio15 Precipitation seasonality, C of V 
Biol6 Precipitation of wettest quarter, mm 
Biol7 Precipitation of driest quarter, mm 
Biol8 Precipitation of warmest quarter, mm 
Biol9 Precipitation of coldest quarter, mm 


matic factors as input data of MaxEnt model to simulate the potential distribution area of the species under 
contemporary and future climate conditions. This study is mainly concerned with the following objects: 
(1) identifying climatically suitable habitats for B. bombina at a home range scale; (2) estimating leading 
climatic requirements (niche) of B. bombina; (3) determining suitable areas for the (re)introduction and/ 
or conservation of the species; (4) forecast the distribution in a future climate, based on a climate change 
projection model for 2050. 

Geographical biases in the occurrence records were dampened by thinning the distribution points with 
OccurrenceThinner 1.03 (Verbruggen et al., 2013). 

Usually, researchers calculate correlation coefficients (e. g. Pearson coefficient) to avoid correlated vari- 
ables and to reduce the effects of multi-colinearity in their models. However, from this type of analysis eco- 
logically relevant variables could be excluded. Burnham and Anderson (1998) have made clear that applying 
correlation analysis in order to find a significant set of predictor variables will most probably expose false cor- 
relations. Fortunately, MaxEnt is able to incorporate complex dependencies between predictor variables, even 
in the presence of correlated variables, non-linearity, bimodality etc. Thus, all covariates were retained for the 
final model. 

The jackknife variable importance feature in MaxEnt was used to assess the relative importance of the 
environmental predictors in the model. We determined the relative importance of variables remaining in the 
model with percent contribution. MaxEnt allows the construction of response curves to illustrate the effect of 
selected variables on probability of occurrence. Identification of the most important predictors, and the analysis 
of the relations between the predictors and predicted habitat suitability, allows the description of the autecology 
of a species. 

We ran the MaxEnt models using the default setting, except for when selecting regularization values. This 
parameter was determined by the application of the small sample corrected variant of Akaike’s Information 
Criterion (AICc) implemented in ENMTools 1.3 (Warren, Seifert, 2011). 

Model performance was assessed using t Warren he average AUC (area under the receiver operating 
curve) score. AUC values >0.9 are considered to have “very good”, > 0.8 “good” and > 0.7 “useful” discrimina- 
tion abilities (Swets, 1988). The logistic output format was used, because it is easily interpretable with logistic 
suitability values ranging from 0 (lowest suitability) to 1 (highest suitability). Better interpretation is made in 
most cases by defining thresholds of habitat suitability. 

With a reference to the classification proposed by Yang et al. (2013), five classes of potential habitats can 
be distinguished: unsuitable habitat (0-0.2); barely suitable habitat (0.2-0.4); suitable habitat (0.4-0.6); highly 
suitable habitat (0.6-0.7); very highly suitable habitat (0.7-1.0). For each model, we calculated the area of the 
optimal distribution, classified as highly or very highly suitable habitat (0.6-1). 
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Results 


In total, our database on B. bombina consisted 2,277 georeferenced point data. After 
reducing the geographical sampling bias in occurrence records their number was thinned 
to 598. Several regularization values were tested: 1, 1.5, 2, 2.5 and 3. The better regulariza- 
tion value (lower value of AICc) was 2. We ran models with 10 bootstrap replicates. The 
resulting potential distribution of B. bombina within its native range under contemporary 
climate (fig. 1) had high AUC values (0.836 + 0.002). In terms of the bioclimatic niche, the 
most suitable for the species areas are located in Poland and Ukraine (fig. 2). In Ukraine 
areas of high or very high suitable habitat potential (> 0.6) occupy around 46 % of the coun- 
try, in Poland this figure is even higher — 65 %. Together these countries accommodate 
87 % of the potential for highly or very highly suitable habitat for B. bombina. 

At the northern border of the distribution of the toad areas of moderate suitability are 
found in southern Latvia. Areas of moderate suitability are found as well in Germany, Den- 
mark and Sweden, but in these countries there are sporadic patches of even highly suitable 
habitat, which should favour the (re)introduction of the species. B. bombina has been suc- 
cessfully reintroduced to Sweden (Andren, Nilson, 1988), but the modelling suggests that 
the species has the chances to spread to a wider portion of the country. 

The MaxEnt model’s internal jackknife test of variable importance showed that the 
environmental variable with highest gain when used in isolation is Bio2 (Mean diurnal 
temperature range), which therefore appears to have the most useful information by itself. 
The environmental variable that decreases the gain the most when it is omitted is Bio3 (Iso- 
thermality), which therefore appears to have the most information that isn’t present in the 
other variables. The percent contribution of Bio2 (29.4 %) too supports the significance of 
this variable for contributing to the MaxEnt model. Another significant (> 10 %) contribu- 
tion to the model is made by Bio8 (Mean temperature of wettest quarter, i. e., predomi- 
nantly summer, 14.9 %). 

Further on, response curves gave an indication of the dependence of the predicted 
probability of presence (i. e., predicted suitability) on each of these variables as well as the 
range under which they reach an optimum of suitability. These probabilities are calculated 
for the range values of one variable, with all other 19 variables set to their average value over 
the set of presence localities. 

Upward trends for variables indicate a positive relationship; downward movements 
represent a negative relationship; and the magnitude of these movements indicates the 





Fig. 1. The potential distribution map for B. bombina under contemporary climatic conditions. The colour 
gradient represents high (red) to low (green) habitat suitability for the species. 
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Fig. 2. The predicted areas of high or very high suitable habitat potential (> 0.6) for B. bombina (A — for con- 
temporary climate conditions; B — for projected climate in 2050). 
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Fig 3. Response curve showing how the logistic prediction changes as the environmental variable Bio2 (Mean 
diurnal temperature range, °C, X-axis) is varied, keeping all other environmental variables at their average 
sample value. The curve shows the mean response of the 10 replicate Maxent runs (red) and and the mean +/- 
one standard deviation (blue). 


strength of the relationship (Baldwin, 2009), essentially they demonstrate biological toler- 
ances. The corresponding response curve for Bio2 shows a clear downward trend (fig. 3), 
representing a negative relationship between increasing values of the Mean diurnal tem- 
perature range and predicted probability of presence (i. e., habitat suitability) for B. bom- 
bina in the locality. 

Under future climate change warmer temperatures are expected to lead to the conver- 
gence of daytime and nighttime temperatures, meaning a reduction of the diurnal tempera- 
ture range (Rohr, Raffe, 2010). These effects of temporal variation in temperature will be 
more pronounced in temperate regions and theoretically may turn out to enhance habitat 
suitability for B. bombina. 

Compared with the area of the most optimal habitat under current climate predic- 
tion, the predictions for 2050 show a clear north and north-west shift of such habitat. It is 
questionable, of course, can the toads track their climate space. This may greatly depend 
on their migratory abilities and how fragmented is the landscape. Once again, in terms of 
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Fig. 4. The potential distribution map for B. bombina under projected 2050 climatic conditions. The colour 
gradient represents high (red) to low (green) habitat suitability for the species. 


the bioclimatic niche, large suitable (> 0.6) for the species areas will continue to persist 
in Poland. Conditions are predicted to considerably improve for the species in Lithuania, 
neighbouring areas of Belarus, as well as in Latvia, where patches of highly suitable habitat 
conditions will appear. Moderate habitat suitability may promote the species to expand as 
far as Estonia and up to the latitude of Saint Petersburg in Russia (fig. 4). 

Optimal habitat is predicted to expand in Sweden, but is likely to shift to the north-east 
from its current location. If so, this may have an impact on measures for reintroducing the 
species. 

In Ukraine, however, areas of potential high or very high habitat suitability are pre- 
dicted to drastically drop from 46 % of the country to 18 %, meaning there could be a net 
extinction of B. bombina at the southern boundary of the home range of the species. In 
the context of these shifts the species’ representation in protected areas can seriously be 
affected, particularly in the south and east of the country. Preemptive measures should be 
undertaken to preserve the toad and its associated habitats in the west and north-west of 
Ukraine, where favourable conditions are predicted to persist. 


Conclusions 


Species distribution modelling has confirmed that vast areas in Poland and Ukraine, in 
terms of the bioclimatic niche, are highly suitable for the Fire-bellied toad, B. bombina. Fur- 
ther to the north areas of moderate suitability are found in a number of countries ranging 
from Denmark to Latvia. The modelling suggests that prospects for the species in Sweden 
could be better than thought. 

Modelling has shown that the Mean diurnal temperature range is an important di- 
mension of the bioclimatic niche of the toad and the expected convergence of daytime and 
nighttime temperatures under global climate change may favour the species. 

Under climate predictions for 2050 there may be a substantial north and north-west 
shift of optimal habitat. Under such a scenario B. bombina populations may suffer mostly 
in the east and south of Ukraine. Optimal habitat is predicted to expand in Sweden, but is 
likely to shift from its current location. This may affect ongoing efforts to reintroduce the 
species to the country. 

Under the modelled scenario the species representation in protected areas throughout 
the home range should be considered, but especially in Ukraine. 
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